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SUMMARY 


The capabilities of a computer program obtained from the Spalding group 
at the Imperial College (London) and further developed by NASA are explored, 
and computed results are compared with data. The comparisons are restricted 
to two-dimensional (2-D) flows; that is, to mixing between parallel streams, 
plane or axisymmetric. Subsonic and supersonic flows, ducted and nonducted, 
reacting and nonreacting, are considered. Evaluation of models used for tur- 
bulence and chemical reaction form an important part of the study. 

Constants in the k - e turbulence model, which produces mixing in good 
agreement with data, are the same for all calculations, but good initial pro- 
files of turbulence kinetic energy k and dissipation rate of turbulence 
kinetic energy e are necessary to obtain satisfactory results. To obtain 
agreement with data in reacting flows, adjustment is necessary over an order- 
of-magnitude range of a constant in the eddy-breakup reaction model. Since 
calculated results in agreement with the data can be obtained only by adjusting 
initial profiles in the k - e model and a parameter in the reaction model, a 
true predictive capability is not present. However, the program is useful for 
interpreting and extrapolating the results of experiments. 

Experimental data were largely drawn from the literature, but new data are 
reported for coaxial injection at matched pressure (1 atm or 101.3 kPa) of a 
cold, Mach 2, hydrogen jet into a hot, Mach 2, vitiated airstream. Profiles of 
pitot pressure and gas composition obtained from water-cooled probes are reported 
and compared with theoretical results. 


INTRODUCTION 

A vital part of the effort (ref. 1) to develop a supersonic combustion ram- 
jet engine (scramjet) is the capability for theoretical prediction of flow prop- 
erties in the combustor, where hydrogen is injected into and burned in a super- 
sonic airstream. For complete generality, the prediction method should be able 
to calculate turbulent mixing and chemical reaction in a combustor for arrays 
of hydrogen injectors at various angles to the airstream. The effects of the 
combustor walls on the flow should be included, and an ability to predict the 
effects of ignition sources and flame holders (in premixed flow) is needed. 

Shock waves and expansions must be accounted for. 

Practical limitations of various kinds make difficult the inclusion of all 
these capabilities in a single program. The present program is restricted to 
mixing of parallel streams, subsonic or supersonic, plane or axisymmetric, react- 
ing or nonreacting, ducted or nonducted. Since it is a two-dimensional program, 
three-dimensional effects such as injection from arrays or at angles to the flow 
are excluded. Interaction with walls is handled by using wall functions (ref. 2), 
a device which transmits the influence of a wall through the wall boundary layer 



without calculating details of the wall boundary-layer flow. Turbulence is 
modeled using the k - e model (ref. 2) and chemical reaction is modeled using 
a modified form of the eddy breakup model (ref. 3). 

Many reports have been published in which experimental mixing data for tur- 
bulent flows have been compared with predictions made by using finite-difference 
programs similar to the one described here. Some of these used the k - e tur- 
bulence model (or one of its "cousins" in which the dissipation rate e is 
replaced with a frequency or a scale length) . The collection of papers in ref- 
erence 4 gives a representative sample of such work. 

Few reports are available which give measured composition profiles in tur- 
bulent flames - especially for supersonic flow. Also, there are no satisfactory 
models available for theoretical calculation of chemical reaction in turbulent 
flows. Some new data on composition in a supersonic flame are published here, 
and emphasis is laid on testing the abilities of the chosen turbulence model 
and the chosen reaction model to predict the properties of turbulent reacting 
flow. 


NOMENCLATURE 

a mass fraction 

Cebu coefficient in eddy breakup model 

Cy , C e i , C e 2 , Cg-| , Cg 2 coefficients in turbulence model 
D nozzle diameter, m 

dj jet diameter, m 

f acceleration parameter in turbulence model 

g mean square fluctuation in mass fraction 

h slot height, m 

h total enthalpy, m 2 s” 2 

i parameter in transport equations 

k turbulence kinetic energy, m 2 s“ 2 

& m mixing length, m 

characteristic length for dissipation of k, m 
M Mach number 

m mole fraction 
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p static pressure, Pa 

p re f reference pressure, Pa 

Ptj2 pitot pressure, Pa 

r radius, m 

upper limit of r for mass flow integral 
T temperature , K 

u flow velocity in principal flow direction, ms - *! 

v flow velocity in transverse direction, ms -1 

w chemical reaction rate, kgm“3s-1 

x distance along flow direction, m 

y distance normal to flow direction, m 

e dissipation rate of turbulence kinetic energy, m 2 s“3 

U viscosity, Nsm -2 

p density, kgm"3 

b Prandtl number 

T shear stress, Nm“ 2 

Subscripts : 
cl center line 

e edge 

f fuel 

g fluctuation 

h enthalpy 

J jet 

k turbulence kinetic energy 

m species 

o initial 
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ox oxidant 

t turbulent 

u velocity in flow direction 

e dissipation rate 

Chemical notation: 

Ar argon 

H 2 hydrogen 

H 2 O water 

N 2 nitrogen 

O 2 oxygen 


DESCRIPTION OF PROGRAM 


The two-dimensional (2-D), parabolic computer program is described in ref- 
erences 5 and 6. It is basically the program the contractor delivered but has 
been tailored to fit the needs of the users. Solutions of parabolic partial 
differential equations for transport of momentum, energy, and mass are accom- 
plished by the finite-difference technique of Patankar and Spalding (ref. 7). 
These equations are 


pu 3u + pv 3u 
3x 3y 


1 _ 

y 1 



i 3u\ 

3y/ 



p u 3h + p y lik 
3x 3y 


_1_ 3 Ait y i 3h\ 

y i 3y\Oh 3y/ 



- + /lit 

a h /3y\2 ) \a k 



( 1 ) 


( 2 ) 


pu 3a + py 3a _ 3 A*t yi 3a\ - # 

3x 3y yi 3y\o m 3yj 


(3) 


where i = 0 for plane flow, and i = 1 for axisymmetric flow. 


Two more equations are solved for the transport of turbulence kinetic energy 
k and of the dissipation rate of turbulence kinetic energy e 

pu IK + pv IK. - J_ 3_/Ht yi 3k\ _ ^ t /3u\ 2 _ pe ( 4 ) 

3x 3y yi 3y\a k 3y ) \3y/ 
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(5) 
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Equations are also solved for the transport of fluctuations in fuel and oxidant. 


pu + pv - 1_ 3_A^ 3ff^ 

3x 3y y i 3y\a g 3y j 


Cg1 (S)(i) 2 ' Cg2 (f) pgf 
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( 6 ) 

(7) 


Turbulence viscosity is computed from 

Mt- 

e 


( 8 ) 


The chemical reaction rate^ is computed from 

w b -CEBuPS 1/2 ^|j (9) 

The values of the constants used were those recommended by authors of the papers 
from which the equations were taken (refs. 3, 5, and 6). These values are as 
follows : 


Cy 

Cel 

Ce2 

C g1 

C g2 

Cebu 



Other a values 

0.09 

1.43 

1.92 

2.80 

2.00 

0.53 

1 .0 

1.3 

0.7 


As explained in reference 2, Cy and C e2 are modified for axisymmetric flow 
as follows : 

Cy = 0.09 - 0.04f 
C e2 = 1.92 - 0.0667f 


where 


Y 

AU 


Y ( 9u cl _ 

9u cl 

\ 

2AU\ 3x 

I 3x 

J 


radial width of mixing region 

axial direction velocity difference across width of mixing region 


^The ability to calculate reaction rates was not included in the program 
as received from the contractor; equilibrium chemistry was used there. 
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Equation (9) for the chemical reaction source term w is the same as the 

one given in reference 3, which was concerned with premixed flow. To use it for 

diffusion flames , the value of g 1 / 2 was taken to be the root mean square of 
the fluctuations in fuel or oxidant, whichever was smaller. This assumption was 

necessary, since w must vanish if either fuel or oxidant is not present . 

Unfortunately, the recommended value Cebu = °“53 was not found to give the 
best fits of calculated results with data. It was necessary to vary its value 
over a range of at least an order of magnitude. This point is discussed later 
in the paper. 

The program provides for selection of two-dimensional or axisymmetric geom- 
etry. Boundaries to the flow may be chosen to be either walls or free bounda- 
ries; a symmetry axis can be specified for axisymmetric flow. A choice can be 
made between H 2 injection or premixed H 2 and air. Chemical reaction can be sup- 
pressed, the gases can react at a finite rate, or they can react completely to 
the extent that they are mixed. A choice can be made between a pressure distri- 
bution which varies only in the flow direction and one in which there are nonuni- 
form transverse pressure profiles. Any desired initial profiles of turbulence 
kinetic energy and its dissipation rate can be specified, or they can be calcu- 
lated in the program by using a mixing length model. Free-stream levels of k 
can be specified independently in the two streams. Ignition can be delayed to 
any point downstream of the initial station. 

Although provision is made for the input of initial profiles of seven vari- 
ables, only profiles of temperature, axial velocity, pressure, and H 2 mass frac- 
tion are required; the others (k, e, and v) are optional. Up to 80 profile 
points may be specified. The initial profile points can be used exactly as 
specified or interpolation is available to produce any desired number of points 
up to 80. Interpolation can be linear or by means of a power-law option. Use 
of the power-law option can concentrate grid points at either boundary or in 
the center of the grid. Linear interpolation can apply over the whole grid, or 
in two sections, where the division between sections is at the H 2 jet boundary. 
Total mass flow of H 2 injected is an input and is used to adjust the initial 
grid in the H 2 region to insure that the calculated mass flow of H 2 equals the 
input value. One other option calculates 1/7 -power-law velocity variations in 
the two parts of the flow. 

In normal operation the transverse grid covers the mixing region between 
the H 2 jet and the co-flowing stream and extends a short distance into the free- 
stream flow. As the mixing region expands, departure from zero of property 
gradients (either velocity or temperature) near a boundary causes the grid to 
expand by entraining free-stream fluid, but the total number of grid points is 
invariant; that is, the grid stretches over a larger region of space to accommo- 
date the additional fluid. If the grid reaches a wall, boundary conditions are 
reset to fix that boundary at the wall and no more fluid is entrained through 
that boundary. A corresponding procedure is used if the grid reaches an axis 
of symmetry. Of course, the calculation can begin with the grid already at the 
boundaries. One other interesting mode of operation provides for calculating 
with specified variation of pressure in the principal flow direction, no entrain- 
ment of fluid , wall boundary conditions, and with variable (calculated) flow 
area. 
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Thermodynamic properties are based on data from the JANAF tables (ref. 8) 
and are made available in the program by means of the least-squares-ad justed 
power-law fits described in reference 9. Temperature is obtained by iteration 
to match the known local value of static enthalpy . Density is calculated from 
the gas law. Computer storage is 111 OOOq and the run time is about 5 minutes 
on a Control Data Corporation 6600 machine to go 30 jet diameters downstream 
from the injection point when 61 transverse grid points are used. 

It is difficult to assess the accuracy of the finite-difference procedure 
used to solve the differential equations. A short qualitative discussion of 
accuracy appears in reference 7 (p. 156), but for practical purposes it was 
found to be more satisfactory to compare with results obtained from other cal- 
culations or with measured data. An important requirement for hydrogen injec- 
tion studies is that the mass flow of hydrogen integrated over the cross section 
at any given station be constant; this hydrogen mass flow is found to vary less 
than 1 percent. 


DESCRIPTION OF TEST CASES 

The several test cases described here were selected to demonstrate the 
ability of the program to make useful predictions of the mixing and reaction of 
hydrogen and oxygen under a variety of conditions. Table I lists the character- 
istics of those selected. 


TABLE I.- CHARACTERISTICS OF TEST CASES 


Test 

case 

Supersonic 

flow 

Geometry 

Static i 

temperature 

Boundary 

Reaction 

Source 


Inner 

Outer 

Inner 

Outer 



1 

Yes 

Axisymmetric 

251 

1495 

Axis 

Free 

Yes 

Present report 

2 

Yes 

Axisymmetric 

276 

1140 

Axis 

Free 

Yes 

Reference 10 

3 

No 

Axisymmetric 

300 

300 

Axis 

Free 

Yes 

Reference 1 1 

4 

Yes 

Axisymmetric 

255 

219 

Axis 

Free 

No 

Reference 12 

5 

Yes 

Plane 

254 

1270 

Wall 

Free 

Yes 

Reference 13 

6 

Yes 

Plane 

294 

294 

Wall 

Wall 

No 

Present report 
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Test Case 1 


For test case 1 (Beach, appendix A), the geometry of the supersonic nozzle 
and H 2 injector is illustrated in sketch (a) : 

D = 0.0653 m 



The hydrogen was piped from high pressure bottles and was expanded through 
a supersonic nozzle (not shown) near the end of the injector pipe. A test gas 
which simulated hot air was obtained by burning H 2 and air, replacing the oxy- 
gen, and then expanding through a supersonic nozzle. Data consisted of pitot- 
pressure profiles and composition profiles. Data at each of the four stations 
were taken in separate runs. Values of x/dj for the composition profiles are 
larger than those for the pitot-pressure profiles because of the geometry of the 
sampling probes. This and other details are discussed in appendix A. The test 
conditions for this case are 



Hydrogen jet 

Free stream 

Mach number , M .... 

2.00 

1.90 

Temperature, T, K . . . 

251 

1495 

Velocity, u, m/s . . . 

2432 

1510 

Pressure, p, MPa . . . 

0.1 

0.1 

Mass fraction: 

ajj2 • • • * • 

a C>2 • 

1.000 

0 

0 

0.241 

a N2 ......... 

0 

0.478 

a H 2 0 ......... 

0 

0.281 
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Test Case 2 


For test case 2 (Cohen and Guile, ref. 10), the geometry of the supersonic 
nozzle and H 2 injector is illustrated in sketch (b) : 


D = 0.100 m 
dj » 0.0200 m 

Injector lip thickness * 0.00054 m 



The hydrogen was piped from high pressure bottles and was expanded through 
a supersonic nozzle (not shown) near the end of the injector pipe. A test gas 
which simulated hot air was obtained by burning H 2 and air, replacing the oxy- 
gen, and then expanding through a supersonic nozzle. Data consisted of pitot- 
pressure profiles, composition profiles, and temperature profiles as given in 
reference 10. The test conditions for this case are 



Hydrogen jet 

Free stream 

Mach number, M . . . . 

1.46 

1.86 

Temperature , T, K . . . 

276 

1140 

Velocity, u, m/s . . . 

1877 

1265 

Pressure, p, MPa . . . 

0.09 

0.09 

Mass fraction: 



a H 2 

1.0 

0 

a 02 

0 

0.260 

a N2 

0 

0.590 

a H 2 0 

0 

0.150 
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Test Case 3 


For test case 3 (Kent and Bilger, ref. 11), the geometry of the air nozzle 
and H 2 injector is illustrated in sketch (c) : 

L = 0.305 m 
d. = 0.00762 m 

Injector lip thickness = 0.00020 m 



The circular H 2 injection nozzle was externally and internally contoured 
to produce a very thin edge at the exit. Turbulence intensity in the square- 
sectioned wind tunnel was claimed to be about 0.2 percent. Data consisted of 
profiles of dynamic pressure, composition, temperature, velocity, and turbulence 
intensity as given in reference 11. The test conditions for this case are 



Hydrogen jet 

Free stream 

Mach number, 

M .... 

0.135 

0.043 

Temperature , 

T, K . . . 

300 

300 

Velocity, u, 

m/s . . . 

178 

15.1 

Pressure, p, 

MPa ... 

0.1 

0.1 

Mass fraction 

9 

« 



a H 2 • • • 

ac>2 0 ’ * 

a N2 ... 


1.0 

0 

0 

0 

0.232 

0.768 

a H 2 0 • • * 


0 

0 
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Test Case 4 


For test case 4 ( Eggers , ref, 12) , the geometry of the supersonic nozzle 
and H2 injector is illustrated in sketch (d) : 

D * 0.152 m 



The hydrogen was piped from high pressure bottles and was passed through 
the nozzle at subsonic speed (M = 0.88) into the supersonic airstream (M = 1.3! 
Data consisted of velocity profiles and composition profiles as given in refer 
ence 12. There was no burning because both the hydrogen and the air were cold 
The test conditions for this case are 







Test Case 5 


For test case 5 ( Burrows and Kurkov, ref. 13) , the geometry of the two 
plane mixing layers is illustrated in sketch (e) : 

Injector lip thickness * 0.0076 m 



h * 0.004 m 


Sketch (e) 

Hydrogen was injected at sonic speed through a slot in the lower wall of 
the test section. A test gas which simulated hot air was obtained by burning 
H 2 and air, replacing the oxygen, and then expanding through a supersonic noz- 
zle. Data consisted of pitot-pressure profiles, composition profiles, and total 
temperature profiles as given in reference 13. The test conditions for this 
case are 



Hydrogen jet 

Free stream 

Mach number, M . . . . 

1.00 

2.44 

Temperature, T, K . . . 

254 

1270 

Velocity, u, m/s . . . 

1216 

1764 

Pressure, p, MPa . . . j 

; 

0.1 

0.1 

Mass fraction: 



a H2 • 

a (>2 

a N 2 

a H20 

1 .000 
0 
0 
0 

0 

0.258 

0.486 

0.256 
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Test Case 6 


For test case 6, the geometry of the supersonic duct is illustrated in 
sketch (f ) : 




-*»Air 


0.02 m 

0.01 m 

|4 — 0.02 m ►+< — 0.02 m — I 

7 7 7 y 7 ^ 7 7 7 

a a 3.366° 

Sketch (f) 


0.0194 m 


Air at Mach 1.5 flows in a two-dimensional duct with walls infinitely dis- 
tant on two sides. The distance between the other two walls decreases linearly 
between two given points along the flow direction; everywhere else the walls are 
parallel. A shock wave with origin at the beginning of the contraction propa- 
gates across the flow and reflects from the opposite wall. In similar fashion, 
a rarefaction with origin at the end of the contraction propagates into the flow. 
This is not an experimental test case; comparison is made with results obtained 
from a shock-fitting program (ref. 14). The test conditions for this case are 



Airflow 

Mach number, M . . . . 

1.50 

Temperature, T, K . . . 

294 

Velocity, u, m/s . . . 

516 

Pressure, p, MPa . . . 

0.1 

Mass fraction: 

a H 2 

a 02 ■••■*■•••••• 

a N 2 • 

a H 2 0 

0 

0.232 

0.768 

0 
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CALCULATED RESULTS 


With one exception the experiments described in the test cases have 
appeared in the literature and are not described here in detail. The data of 
test case 1 , in which hydrogen was injected into a co- flowing, supersonic, axi- 
symmetric, heated airstream, were taken by using the equipment and procedures 
described in appendix A. The option for computation of transverse pressure 
profiles was used only for calculating test case 6, which was chosen especially 
to demonstrate this capability. (There are no experimental data for test 
case 6.) 

The purposes of comparing experimental and calculated results were to 
validate the use of the program under a variety of conditions and to determine, 
through experience, any special sensitivities or limitations to its use. The 
ability of the turbulence model to provide eddy viscosities leading to good pre- 
diction of mixing under a wide variety of conditions was of especial interest - 
as was also the ability to account for the effects of chemical reaction. 


Test Case 1 

Because test case 1 (Beach, appendix A) to a substantial degree simulates 
a portion of the reacting flow in a hydrogen- fueled supersonic combustor, and 
because the measured data are previously unpublished, the analysis is more com- 
plete than those for the other cases. Initial profiles at x/dj = 0.33 are 
given in figure 1 for temperature, pitot pressure, velocity, composition, turbu- 
lence kinetic energy, and its rate of dissipation. Static pressure was assumed 
to be 1 atmosphere throughout the flow field; the composition consisted of hydro- 
gen in the jet and vitiated air outside. The initial static temperature and 
velocity profiles were based on the known total temperatures in the jet and 
free stream and on the measured profile of pitot pressure. 

The flow calculations began at x/dj s 0.33, the point at which measure- 
ments closest to the injector were made. It should be noted that the interior 
and exterior diameters of the injector were 0.0065 and 0.0095 m, respectively; 
the lip thickness was 0.0015 m. A small recirculation region probably was pres- 
ent near the lip. If any shock waves were generated by interaction between the 
jet and the surrounding flow, their effects on the calculated flow properties 
are not accounted for. Since the jet and free-stream static pressures were 
carefully matched in the experiments, it is believed that no significant effects 
caused by such shock waves were present. 

Axial variations of pitot pressure and hydrogen mass fraction are shown 
in figure 2. The reasonable agreement between theory and experiment is a good 
^indication that mixing is calculated correctly. Since mixing is produced by 
turbulent diffusion in the direction transverse to the axis, and since the mag- 
nitude of the diffusion term is proportional to the turbulent viscosity modeled 
by equation (8), it is important to have a good procedure for determining the 
initial profiles of k and its dissipation rate e. Experimental profiles of 
k and e were not available; therefore, Prandtl 1 s mixing length theory was 
used in the manner described in appendix B. Unfortunately, adjustment of the 
dissipation length parameter Jl e is necessary to achieve good agreement with 
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the data. To illustrate this point , curves are shown in figure 2 for three 
values of 


Since the location of streamlines in a reacting flow can depend strongly 
on the amount of fuel that reacted, it is also important to examine the radial 
variation of the integrand in the mass conservation relation : 


Mass flow = 



pur dr 


( 10 ) 


In figure 3 experimental 2 * and theoretical values of pur at four axial stations 
are compared. The complete-reaction curve corresponds to an infinitely fast 
reaction; that is, all fuel and oxygen that mix are completely burned. The 
quantity Cebu is the constant that appears in equation (9); thus, the curves 
identified by a value of Cebu were calculated by using equation (9). The 
no-reaction curve was obtained by using frozen chemical composition. Best 
agreement is obtained by using Cebu = 0.10; the other curves appear on the 
plots to give the reader a feeling for the range of values to which the calcu- 
lations are sensitive and to indicate roughly the degree to which reaction is 
complete. More precisely, for this case (Cebu = 0.10), the computed reaction 
efficiency (defined as the ratio of hydrogen burned to hydrogen macroscopically 
mixed) is 0.20 at x/dj = 6.56, 0.28 at x/d< = 13.8, 0.30 at x/dj = 20.0, and 
0.32 at x/dj = 26.2. 


It will become apparent as other test cases are discussed that Cebu is 
not a constant; its value must be determined by analysis of measured data like 
that shown in figure 3. Profiles of pitot pressure and composition calculated 
using Cebu = are Shown in figure 4. No-reaction and complete-reaction 

curves are shown for comparison. The data points for O 2 are the most sensitive 
indicator of the degree to which reaction is complete, and these show reasonable 
agreement with the Cebu = 0 • 1 0 curves . The data points for H 2 O do not agree 
as well and consistently indicate more water formed than was predicted by the 
Cebu = 0.1° calculation. However, the pitot-pressure data do not reinforce 
this observation; they indicate less reaction than the Cebu = 0- 10 calculation. 
A plausible explanation for this result is that additional reaction takes place 
in the gas sampling probe. As discussed in reference 15, if this reaction 
occurs, the gas sample composition is greatly altered, but the pitot pressure 
is not substantially affected (3 percent to 12 percent maximum). 


It should be noted that the eddy breakup model assumes one global reaction 
(O 2 + 2 H 2 ■*► 2 H 2 O) and makes the rate proportional to the rate at which large 
turbulent structures break down into small ones. A model based on a system of 
chemical kinetic equations would be able to produce H 2 O faster in the center 
of a flame where temperatures and concentrations are high and would produce it 
slower at the edges. Such an ability might produce calculated curves which 
agree better with the data. 


2 Density and velocity profiles were calculated from the experimental data 

on pitot pressure and composition. 
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It should also be remembered that in turbulent flames like the one being 
examined, unmixedness can play a major role. Since the eddy breakup model is 
based on the concept of unmixedness, it is likely to be superior in some situa- 
tions to chemical kinetic models which do not account for unmixedness. It would 
be desirable to have a model which can account for both chemistry and unmixed- 
ness, but such a model is beyond the scope of this report. 


Test Case 2 

Calculations were made for test case 2 (Cohen aftd Guile, ref. 10) because 
the test conditions were similar to those of test case 1 . In both experiments 
hydrogen was axially injected at supersonic speed into a supersonic stream of 
vitiated air. Geometry was similar and composition profiles were measured down- 
stream in both cases. The most evident differences were in the temperatures of 
the airstreams (1140 K for this case and 1495 K for test case 1) and in the 
sizes of the jets (dj = 0.0200 and 0.009525 m) . 

The plots in figures 5 to 7 should be compared with those in figures 2 to 
4. In both cases the initial value of was adjusted to give reasonable 

agreement with the measured data. The distributions of mass flow (figs. 3 and 
6) are similar and show that the amount of reaction is small for both cases, 
although more reaction appears to have occurred in the hotter flow (case 1). 

For the case 2 finite-rate calculation using Cggy = 0.05, the computed reaction 
efficiency is: 0.14 at x/dj = 5.1, 0.15 at x/dj = 8.9, 0.16 at x/dj = 12.7, 
and 0.16 at x/dj = 17.8. 

The pitot-pressure data in figure 7 and the mass flow distributions in 
figure 6 agree best with the no-reaction curves, but the H 2 O data and the O 2 
data indicate that some reaction occurred. The value Cgg^ =0.05 is used to 
demonstrate the effect of a small amount of reaction and calculated results 
obtained with it are not offered as a best fit to the data. 


Test Case 3 

The experiment for test case 3 (Kent and Bilger, ref. 11) was again axial 
injection of hydrogen into a co- flowing stream of air. Unlike the two previous 
cases, the Mach number was low in both streams. Also, the airstream was cold 
and did not contain the large amount of water present in vitiated air. Great 
care was taken in the experiment to contour both the hydrogen and the air noz- 
zles so as to produce low initial turbulence; turbulence intensity was claimed 
to be about 0.2 percent in both the jet and the free stream. 

Figure 8 shows the axial variation of several properties for case 3. 

Except for the momentum flux at the outer edge of the jet p e u e 2/(p e u e 2) 0 , all 
the properties shown are on the center line. The data extend much farther down- 
stream than in cases 1 and 2 (x/dj = 1 60 as compared with x/dj = 30). Agree- 
ment between theory and experiment is good. The largest discrepancies are in 
the temperature, and even here the trends are correct. 
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The mass-flow distributions shown in figure 9 are strikingly different 
for the no-reaction and the complete-reaction assumptions. Unlike the two pre- 
vious cases, test case 3 data agree with the assumption of complete reaction 
rather than that of little or no reaction. The disagreement between theory and 
data for r t 0.05 and x/dj ^ 80 is thought to be caused by the inability 
of the k - £ turbulence model to predict the correct values of turbulent 
viscosity in the outer region of the flow where velocity gradients are small. 
Support for this conclusion is contained in the results of an experimental inves- 
tigation into how large the errors are in the prediction of turbulent viscosity 
for various flow geometries and in different parts of a flow (ref. 16). The 
results of the investigation show that for axisymmetric geometry (and plane 
geometry also), the viscosity is predicted well in the central part of the mix- 
ing layer, but not in the low shear regions at the edges. This kind of flaw in 
the model can be tolerated because the most important viscosity effects occur 
in the central part of the mixing layer. 

The radial profiles shown in figure 10 demonstrate again that the theoreti- 
cal and experimental results agree well for this case. The agreement is better 
at x/dj = HO than at x/dj = 120 and is better for small values of r/dj than 
in the outer part of the flow. Best results were obtained with Cggy = 0.53. 
Reaction is nearly complete, as figure 9 also indicated. The computed reaction 
efficiency is about 0.90 at both x/dj = 40 and at x/dj = 120, The properties 
of this low-speed diffusion flame are predicted well except in the outer flow 
region at large distances from the injection point. 


Test Case 4 

This experiment (test case 4 (Eggers, ref. 12)) is another in which H 2 was 
injected into a supersonic airstream. There was no reaction, for both gases 
were cold and there was no ignition. Figure 11 shows the decay of center-line 
velocity as a function of distance from the injection point and also shows the 
decay of the H 2 mass fraction along the center line. 

The profiles of mass flow are shown in figure 12 for four stations. Both 
no-reaction and complete-reaction curves are shown, even though case 4 was a non- 
reacting experiment. This is done to demonstrate a similarity between cases 3 
and 4 as compared with cases 1 and 2. The calculated shift of radial position 
between the no-reaction and complete-reaction curves is much larger for cases 3 
and 4. This difference is due to the large density mismatch between the H 2 and 
air layers for cases 3 and 4. 

Consider the effect of density and velocity profiles on the pur plots. 

If Pu values were constant across the mixing layer, the plot would be a 
straight line. In the initial stages, before much mixing occurs, the curve 
consists essentially of two straight lines which correspond to the two free- 
stream conditions. The following tables are helpful for comparing the pur 
plots of the four cases: 
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Case 

Air layer 

P 

u 

pu 

1 

0.203 

1510 

307 

2 

.257 

1265 

324 

3 

1.17 

15 

18 

4 

1.61 

394 

634 


Case 

Ratios (Air/H 2 ) 

pu 

P 

1 

1.29 

2.4 

2 

2,16 

3.2 

3 

1.21 

14.3 

4 

6.16 

16.7 


Case 

h 2 

layer 

P 

u 

Pu 

1 

0.098 

CM 

OO 

■=r 

CM 

238 

2 

.080 

1877 

150 

3 

.082 

178 

15 

4 

.096 

1074 

103 


The air/H 2 ratios of Pu for cases 1 and 3 are both near unity. Thus, 
the slopes of the two parts of the no-reaction curves should be nearly the same. 
The ratios of the slopes are about 2 for case 2 and about 6 for case 4. Crude 
confirmation of these relationships can be had by examining figures 3, 6, 9, 
and 12. 

As was remarked earlier, the density ratios are larger for cases 3 and 4 
than for cases 1 and 2; the outer layers are cold air instead of hot air. When 
reaction occurs, density drops in that part of the cold air near the flame front 
and the streamlines expand to satisfy mass conservation. The effect would be 
the same if the H 2 stream and air stream were interchanged; the radial shifts 
would still be away from the center line and they would be larger for larger 
density ratios . 

It is apparent that Pur plots can be used to determine the amount of 
reaction present, even with quite crude values of measured composition. This is 
a very useful capability, since it is often difficult to be sure that the mea- 
sured sample profiles are representative of the flow being sampled. 

The center-line data for case 4 have been compared previously (ref. 17) 
with calculations made by using the same k - e turbulence model used here. 

The agreement was better than that in figure 11, but profiles were not reported. 
In the present calculation, close attention was given to obtaining good agree- 
ment with the profiles, as shown in figures 12 and 13, but at the cost of some- 
what poorer agreement with center-line data. The agreement is better for veloc- 
ity profiles than for the H 2 mass fraction profiles, but for both, the results 
are acceptable, and the trends indicated are correct. 


Test Case 5 

Test case 5 (Burrows and Kurkov , ref. 13) is the only mixing experiment 
considered in which H 2 was not injected along the center line of an axisymmetric 
airstream; injection in this case was from a slot in a step on the wall of a rec- 
tangular duct. As shown in figure 14, agreement between theory and experiment 
is good when a nearly complete reaction is assumed. In fact, the only place in 
the figure where the finite rate and complete reaction curves differ appreciably 
is near the peak of the curve for mole fraction of water. 
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The agreement is considerably better than that obtained by the authors of 
reference 13. Their calculated profile shapes , including the peak of the water 
profile, matched those of the data well , but were misplaced in the cross-stream 
direction. They discussed the discrepancy themselves and suggested two possible 
explanations: (1 ) that the use of equilibrium, rather than finite-rate chemis- 

try, was responsible for the difference; and (2) that effective viscosity was 
higher in the flame region than their model predicted. They obtained excellent 
agreement between theory and experiment when they examined data obtained in a 
nonreacting experiment similar to the reacting one. 


Test Case 6 

Test case 6 was chosen to illustrate the program's capability for calculat- 
ing transverse pressure profiles. No suitable measured data were available, but 
a good test case was obtained by analyzing the chosen flow field by using the 
shock fitting program described in reference 14. 

Air at room temperature and Mach number 1.5 flows between two plates. All 
property profiles are initially uniform. Between the two given points, the dis- 
tance between the plates decreases linearly; elsewhere, the distance between 
them is a constant. A shock wave, with origin at the first point, moves into 
the flow and is reflected from the opposite wall. In similar fashion, a rare- 
faction moves into the flow from the second point. 

The solid lines in figure 15 show the shock fitting results and the other 
lines show the results obtained with the parabolic marching program. Obviously, 
steep pressure gradients, like the shock waves in figure 15(a), are poorly rep- 
resented. However, the general shape and magnitude of the pressure profiles 
are obtained, and both the shock wave and the rarefaction region are easily 
identified. 

It is interesting to note how the diffuse shock fronts calculated by the 
program cause premature pressure changes at the walls. (See fig. 15(b).) No 
premature pressure change is calculated on the lower wall at y = 0.01 m or at 
y = 0.03 m because these changes were initiated by sudden changes of slope in 
the lower wall. In the program, flow area is caused to match duct area by 
adjusting the average pressure in a one-dimensional approximate treatment. This 
is a better approximation in subsonic flow than in supersonic flow, because the 
finite time required for pressure changes to propagate is ignored in the one- 
dimensional treatment. For an example of this effect, see figure 15(b), in 
which pressure at a wall should be constant except when a shock wave or a rare- 
faction passes. The overall raising of the pressure which can be seen in the 
calculated results is caused by accounting in the one-dimensional treatment for 
boundary-layer growth at the walls. Because of the tendency of shock and expan- 
sion waves to be diffused by mixing and reaction, detailed calculation of their 
propagation in supersonic combustors may not be warranted for many applications. 
The extent to which transverse pressure gradients can be accounted for in the 
present program should be useful in such cases. 
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CONCLUDING REMARKS 


The two-dimensional computer program described here is a useful tool for 
analysis of a variety of turbulent flows - namely, mixing between parallel 
streams, plane or axisymmetric , subsonic or supersonic, ducted or nonducted, 
reacting or nonreacting. The test cases demonstrate that this program can be 
used to calculate properties of such flows with reasonable accuracy. They also 
demonstrate that good modeling of both turbulence and reaction effects is neces- 
sary for calculation of flow properties in turbulent flames. The empirical con- 
stants in the turbulence model were the same for all calculations (assumed uni- 
versality of the model ) , but the initial profiles of turbulence kinetic energy 
k and dissipation rate of turbulence kinetic energy e were found to strongly 
influence the results. The eddy breakup model (modified) for calculating the 
rate of chemical reaction produced composition profiles which in most cases 
closely resembled the measured profiles, but was able to do so only by adjusting 
an empirical constant to cause calculated profiles of mass flow distribution to 
match those measured. Adjustment of this constant controlled the amount of 
heat released by chemical reaction and thereby determined the spatial location 
of streamlines in the flow through the change in gas density profiles. 

The two-equation turbulence model appears to be generally adequate for 
modeling turbulent viscosity in axisymmetric and two-dimensional plane flows, 
except in regions of low shear. Plots of mass flow distribution are very useful 
for deciding how much reaction occurred in a given experiment. The eddy breakup 
model (modified), which is valid only if a flame is diffusion controlled, gave 
reasonable results for the cases reported, but should be replaced by a chemical 
kinetic model, if one can be developed, which adequately accounts for the effects 
of turbulence. 


Langley Research Center 

National Aeronautics and Space Administration 
Hampton, VA 23665 
February 17, 1978 
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APPENDIX A 


DATA FOR TEST CASE 1 

In order to provide additional data for the evaluation of analytical models 
under development, a mixing-reacting experiment was conducted at the Langley 
combustion test stand. The experiment consisted of the injection of a super- 
sonic jet of hydrogen coaxially into a high-temperature , vitiated airstream. 


Apparatus and Procedure 

The test gas for the experiment was generated by burning hydrogen in air 
and replenishing the volumetric oxygen content. The heater which accomplishes 
this is described in detail in reference 18; it is Capable of producing stagna- 
tion temperatures to 2800 K and stagnation pressures to 3.1 MPa. Extensive 
facility calibration has indicated virtually 100 percent combustion efficiency 
for the heater, and the primary test gas constituents are therefore oxygen, 
nitrogen, and water vapor. 

The test gas nozzle, which is constructed of stainless steel, is cooled 
by approximately 6 kg/sec of water. The exit diameter is 6.57 cm, and the exit 
Mach number is approximately 1.9. A stagnation pressure of 0.067 MPa was 
required to achieve the desired 1 -atmosphere static-pressure condition at the 
nozzle exit. Stagnation temperature was nominally 2250 K in order to simulate 
flight enthalpy environments in the Mach 7 regime. The resulting static temper- 
ature and test-gas flow rate were approximately 1500 K and 0.980 kg/sec. Nomi- 
nal test gas constituents were oxygen (20?), water vapor ( 38 ?), and nitrogen 
(42?) by volume. 

The hydrogen injector at the nozzle center line is a 0.95-cm-diameter 
stainless-steel tube with a nozzle insert that produces a 5° exit half-angle. 

The exit Mach number is nominally 2, and the lip thickness approximately 1.5 mm. 
In these tests cooling was provided by the injectant which originated from a 
supply at ambient temperature. Hydrogen stagnation pressure was adjusted to 
provide a matched-pressure condition at the injector and test gas nozzle exits. 
Estimated temperature rise for the hydrogen passing through the injector was 
150 K; this gave a hydrogen stagnation temperature of 450 K, and a hydrogen 
mass flow rate of 0.0084 kg/sec. 

Primary in-stream measurements for the experiment were pitot pressure and 
composition profiles. Details of the pitot probe design are given in refer- 
ence 19. Briefly, the probe tip is a 20° half-angle cone with a 0.79-mm ori- 
fice. The maximum diameter at the base of the probe tip is 0.92 cm. A contin- 
uous traverse of the flow was made at a rate of 0.5 cm/sec. Pressure surveys 
were taken at x/dj locations *of 0.33, 6.56, 13.8, 20, and 26.2; the data 
appear in figures 1 to 4 and in table At. Gas samples were taken with a wedge 
gas sampling probe discussed in reference 20. It has a 20° included angle with 
the sampling surface mounted at 2° to the flow; ports 0.5 mm in diameter are 
drilled 2.4 mm apart. The wedge leading edge has a 0.38-mm radius to reduce 
stagnation-point heating. Gas samples were acquired at x/dj locations of 
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TABLE A1 PITOT-PRESSURE SURVEYS 
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APPENDIX A 


8.26, 15.5, 21.7, and 27. 9, and the results are shown in figures 1 to 4 and in 
table A2 . The difference in axial location between the pitot and gas sample 
surveys is caused by the fact that the wedge-probe and pitot-probe leading edges 
were at the same location, but the sampling ports are 1.7 cm behind the wedge 
leading edge. 

Gas samples were collected in 7 5 -cm3 cylinders (with valves at each end) 
by evacuating them, purging with the sample gases during a run, and finally 
filling them to an appropriate pressure level with the sample gases. Total 
sampling time for a set of nine samples was 8 seconds. Analysis was by gas 
chromatograph, and only dry samples of nitrogen, hydrogen, oxygen, and helium 
were examined; 3 helium was used to trace the replenishment oxygen in order to 
differentiate it from the oxygen in air. Water concentrations were deduced 
from the dry samples by a data-reduction technique described in reference 19. 

A quantitative indication of the validity of this approach is given by compar- 
ison of the reduced free-stream levels of nitrogen, oxygen, and water vapor as 
compared with the bulk values in the heater. (See fig. 4.) 

It should be noted that facility and instrumentation limitations allowed 
only one profile to be obtained in any given burner firing. The nine surveys 
(five pitot, four composition) were therefore obtained in nine different runs, 
and some differences in test gas and hydrogen flow rates were inevitable. How- 
ever , these differences were slight , with maximum deviations of 1 . 1 percent , 

3.8 percent, and 2.2 percent from previously mentioned heater stagnation pres- 
sure, stagnation temperature, and hydrogen fuel flow rate, respectively. The 
stagnation temperatures were not measured directly but were computed from heater 
performance data and the measured flow rates of hydrogen, oxygen, and air to 
the heater. 


Data Integrity 

A traditional indicator of data integrity for mixing flow fields has been 
the comparison of mass flows calculated from measured composition profiles with 
flow rates measured during the experiment. For these free-jet data, the only 
meaningful comparison is for injected hydrogen. Integrations utilizing the 
composition, pitot pressure, and an assumed uniform static pressure of 1 atmo- 
sphere were made at the four sample locations. It was assumed that the closest 
pitot profile to each wedge location was valid at the wedge location. Errors 
resulting from the integration process ranged from 25 to 33 percent. These are 
large errors, although magnitudes of 20 percent are typical for this type of 
flow field. There are several potential sources for the discrepancies. 

Data accuracy of pitot pressure measurements in flows in local thermody- 
namic equilibrium is typically very reliable, and results from the dry samples 
analyzed by the chromatograph are good to within ±1 percent. However, it is 
possible that the probes complete an on-going reaction, and that the constitu- 
ents reaching the sample bottle are not those which were present ahead of the 
probe. If account is taken for this effect by the unreacting of part of the 


3The water was removed before analysis because condensation in the sampling 
tubes would invalidate any attempt to include it in the analysis. 
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TABLE A2 , - SPECIE MASS FRACTIONS 


(a) x/dj = 8.26 


(b) x/dj = 15.5 


m 

Mass fraction for - 

H 2 0 

h 2 

02 

n 2 

Ar 

-0.275 

0.232 

ma 

0 


0.003 

.025 

.141 

.742 

0 

.115 

.002 

.225 

.195 

.642 

0 

.160 

.003 

.475 

.466 

.094 

.025 

.408 

.007 

.725 

.300 

0 

.229 

.462 

.008 

.975 

.267 

0 

.259 

.465 

.008 

1.225 

.257 

0 

.268 

.466 

.008 

1.475 

.247 

0 

.276 

.468 

.008 

1.725 

.244 

0 

.278 

.469 

.008 


r/dj 

Mass fraction for - 


H 2 0 

h 2 

0 2 

n 2 

Ar 

-0.283 

0.315 

0.423 

0 

0.257 

0.005 

-.033 

.215 

.604 

0 

.177 

.004 

.227 

.263 

.517 

0 

.216 

.004 

.477 

.450 

.163 

.004 

.376 

.007 

.727 

.393 

.008 

.137 

.453 

.008 

.977 

.332 

0 

.198 

.461 

.008 

1 .227 

.322 

0 

.207 

.463 

.008 

1.477 

.274 

0 

.250 

.467 

.008 

1.727 

.270 

0 

.253 

.468 

.008 


(c) x/dj = 21.7 


■ 

Mass fraction for - 

H 2 0 

h 2 

0 2 


Ar 

-0.265 

0.353 

0.353 

0 

0.288 

0.005 

-.015 

.312 

.428 

0 

.255 

.005 

.235 

.341 

.374 

0 

.280 

.005 

.485 

.431 

.205 

0 

.357 

.006 

.735 

.440 

.037 

.078 

.437 

.008 

.985 

.362 

.002 

.169 

.459 

.008 

1.23 

.311 

0 

.217 

.463 

.008 

1.49 

.274 

0 

.250 

.466 

.008 

1.735! 

.269 

0 

.255j 

.468 

.008 


(d) x/dj =27.9 


r/dj 

Mass fraction for - 


H 2 0 

h 2 

°2 

n 2 

m 

-0.267 

0.403 

0.262 

0 

0.329 


-.017 

.369 

.323 

0 

.303 

.005 

.233 

.394 

.275 

0 

.324 

.006 

.483 

.429 

.162 

.024 

.377 

.007 

.733 

.445 

.527 

.065 

.429 

.008 

.983 

.380 

.006 

.149 

.456 

.008 

1.233 

.347 

.001 

.183 

.461 

.008 

1.483 

.301 

0 

.224 

.466 

.008 

1.733 

.263! 

0 

.260 

.468 

.008 
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fuel, a 5-percent increase in hydrogen mass flow results from 50-percent local 
unreactedness and 13 percent from 90-percent unreactedness. 

Another possible error source is probe alignment. If the assumption is 
made that the wedge probe was on a chord displaced 0.16 cm off the diameter, 
increases in computed hydrogen mass flows range from 15 percent at x/dj = 8.26 
to 7 percent at x/dj = 27.9. It is difficult to assess the errors produced by 
sampling pitot pressures and composition at slightly different locations, but 
they are believed to be small. 

The assumption of uniform static-pressure profiles is obviously not correct 
even though the pressures of hydrogen and test gas are initially matched. The 
injector has a conical nozzle, and the lip thickness of the injector is large 
enough to cause aerodynamic disturbances. Nevertheless, if a uniform increase 
of 10 to 15 percent is made in the static pressure, the computed mass flow 
increases only 2 to 3 percent. 


Other possible sources of error include the heat loss in the heater and 
nozzle (change in mass flow is 2 percent maximum), and the fact that oxygen 
and helium may not diffuse together (less than 2 percent error). 

Unfortunately, hardware problems eliminated the possibility of examining 
and correcting these potential error sources. In spite of this, the data are 
believed to be representative and meaningful and are appropriate for the eval- 
uations being made in this report. 
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PROCEDURE FOR DETERMINATION OF INITIAL PROFILES OF k AND e 

Turbulent viscosity is modeled by using equation (8) in combination with 
equations (4) and (5). The initial profiles of k and £ used to begin cal- 
culations with these equations have a marked influence on the amount of mixing 
calculated. Since measured profiles of k and e are seldom available, the 
following procedure has been devised for use in their absence. 

First, an estimate is made of background levels of turbulence in each 
stream; these are put in as constant values of k to which shear-generated 
peaks are added. The shear-generated peaks at the initial station are calcu- 
lated by using a mixing-length equation; the result of summing the two contri- 
butions is 

k = ^0^3 + ^background 

where is the dissipation length scale, 

E cb2) 

=• e 


The length scale is chosen equal to 10 times an appropriate physical dimen- 

sion such as the thickness of the splitter plate in a two-dimensional free shear 
layer experiment; it is then adjusted to obtain best agreement with observed 
quantities downstream. Given the k profile from equation (B1), the e pro- 
file is calculated from equation (B2). 

The viscosity p, calculated from equation (8), is a turbulent viscosity. 

In the program the laminar viscosity is computed and is added to p, so that 
the viscosity used in all transport equations is the sum of the laminar and 
turbulent viscosities. Usually, the laminar viscosity makes a negligible con- 
tribution to the total. 

Equation (B2) is a definition. Equation (B1) is derived from Prandtl's 
mixing length equation as follows: 

p = pS, m 2/|u\ (B3) 


This equation is appropriate for flows in which the production term for k is 
equal to the dissipation term, and experiments have shown that for such flows 
(ref. 21), 



0 . 3pk 


(B4) 
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By combining equations (8) and (B2) , 
y = pk 1/,2 Jl e 

From equations (B3) to (B5), 


(B5) 

(B6) 


k = (B7) 

\0.3 3y/ 

Equations (B1) and (B2) are adequate to define initial profiles of k 
and e based on the Prandtl mixing length concept; the only information needed 
is a characteristic dissipation length Jt e and a velocity profile. Experience 
with this approach has been good, provided some experimental knowledge is avail- 
able about the downstream flow so that the value of & e can be adjusted. A 
more satisfactory procedure would be to find il e by referring to a correlation 
plot in terms of known physical quantities so that flow-field predictions could 
be made reliably without reference to measured downstream properties. Unfortun- 
ately, no such relation has been found. 
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Figure 2.- Axial variation of pitot pressure and hydrogen mass fraction on the 
center line for case 1. Cggu =0.10. (H 2 in jet; vitiated air outside.) 





Figure 3*~ Distribution of mass flow for case 1. (H 2 in jet; vitiated 

air outside.) 
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Figure 4.- Profiles of composition and pitot pressure for case 

vitiated air outside.) 
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Figure 4.- Continued. 
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Figure 6.- Distribution of mass flow for case 2. (% in jet; vitiated 

air outside.) 
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Figure 7.- Profiles of composition and pitot pressure for case 2. (H 2 in jet 

vitiated air outside . ) 
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Figure 10.- Profile of temperature, composition, and dynamic pressure 
for case 3. (H2 in jet; air outside.) 
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Figure 10.- Concluded. 
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Figure 11.- Variation of properties on the center line for case 4. (H 2 in jet 

air outside. There was no chemical reaction in this experiment.) 
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(b) x/dj = 9.58. 


Figure 13.- Continued. 














Figure 14,- Profiles of composition, pitot pressure, and mass flow for case 5 
at x/h =89. (H 2 in jet; vitiated air outside.) 






(a) Profiles at x = 0,243 m, x = 0.0356 m, and x s 0.0480 m 
Figure 15.- Pressure distributions for test case 6. 
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